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MULTI-FLUID  PLASMA  MODEL 

AFOSR  Grant  No.  FA9550-11-1-0167 


U.  Shumlak 

Department  of  Aeronautics  and  Astronautics 
Aerospace  &  Energetics  Research  Program 
University  of  Washington 


Abstract 

A  physics-based  algorithm  is  developed  based  on  the  multi-fluid  plasma  model  de¬ 
rived  from  moments  of  the  Boltzmann  equation.  The  model  includes  evolution  equa¬ 
tions  for  the  electromagnetic  fields,  electron  fluid,  ion  fluid,  neutral  fluid,  and  any 
additional  species.  The  large  mass  difference  between  electrons  and  ions  introduces 
disparate  time  and  spatial  scales  and  requires  a  numerical  algorithm  with  sufficient 
accuracy  to  capture  the  multiple  scales.  In  addition,  the  characteristic  time  scales  for 
the  electromagnetic  fields  is  much  shorter  than  the  time  scales  of  the  ion  and  neu¬ 
tral  fluids.  The  physics-based  computational  algorithm  solves  fluid  models  for  each 
plasma  species  that  are  appropriate  for  the  expected  physical  behavior,  by  combining 
5 A- moment  and  13 A- moment  fluid  models  for  multicomponent  plasmas.  The  numer¬ 
ical  discretization  is  developed  specifically  to  capture  the  expected  physical  behavior 
by  combining  high-order  continuous  and  discontinuous  spatial  representations  of  the 
solution  and  implicit  time-advance  methods  to  accurately  capture  the  fast  and  slow 
dynamics.  The  physics-based  computational  algorithm  has  also  been  extended  to  solv¬ 
ing  continuum  kinetic  plasma  models.  Solving  Maxwell’s  equations  has  been  improved 
by  using  a  parabolic  modification  to  ensure  the  errors  of  divergence  constraint  equa¬ 
tions  are  properly  removed  and  handled.  Nonreflecting  boundary  conditions  using  a 
lacunae-based  method  have  been  implemented  and  provide  higher  solution  fidelity  for 
open  boundary  problems. 

1  Project  Description 

Plasmas  are  essential  to  many  existing  and  emerging  technologies  that  are  important 
to  the  Air  Force,  industry,  and  general  science.  These  applications  include  high  power 
microwave  devices,  plasma  actuation  of  airstreams,  drag  reduction  for  hypersonic  vehi¬ 
cles,  advanced  space  propulsion,  weapons  effects  simulations,  radiation  production  from 
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fusion,  etc.  Advancing  plasma  technologies  requires  a  fundamental  and  accurate  under¬ 
standing  of  the  underlying  physical  effects  and  the  resulting  plasma  dynamics.  Plasma 
dynamics  inherently  involve  complex  physical  phenomena  because  the  dynamics  are 
affected  by  short-range  (collisions)  and  long-range  (electromagnetic  fields)  forces.  The 
great  utility  of  plasmas  stems  from  these  multi-scale  forces.  In  general,  plasmas  fall 
into  a  density  regime  where  they  exhibit  both  collective  (fluid)  behavior  and  individual 
(particle)  behavior.  The  intermediate  regime  complicates  the  analytical  and  compu¬ 
tational  modeling  of  plasmas.  Accurate  computational  modeling  of  plasmas  requires 
innovative  numerical  algorithms  that  are  tuned  to  solve  the  specific  physical  models 
that  capture  the  multi-scale  and  multi-physics  phenomena  that  are  present  in  plasmas. 


1.1  Plasma  Models:  MHD,  Kinetic,  PIC 

Understanding  and  predictability  of  plasma  behavior  has  been  significantly  advanced 
through  the  development  of  reduced  plasma  models  and  their  numerical  solution.  The 
most  common  reduced  plasma  model  is  the  magnetohydrodynamic  (MHD)  model, 
which  describes  the  plasma  as  a  single-fluid.  While  the  MHD  model  has  been  success¬ 
ful  in  many  applications,  [1-6]  more  complex  effects  require  more  complete  physical 
models. 

The  most  complete  continuum  model  for  plasma  is  described  using  kinetic  theory 
where  each  species  a  of  a  plasma  is  described  by  a  time-dependent  distribution  function 
/a(x,v,  t)  in  physical  and  velocity  space.  The  evolution  of  the  distribution  functions 
is  described  by  the  Boltzmann  equation 


dfa 

dt 


+  V 


OX 


v  x  B) 


dfa 

<9v 


dfa 

dt 


(1) 


The  plasma  is  composed  of  ion  and  electron  species  and  possibly  additional  species  for 
neutrals  or  impurity  ions. 

The  collision  term  on  the  right-hand  side  of  the  Boltzmann  equation  accounts  for 
changes  to  fa  due  to  short-range  interactions  within  species  a,  between  species  a  and 
species  /3,  and  between  species  /3  and  species  7.  The  three  collisional  interactions  refer 
to  thermodynamic  equilibration  within  a  species,  thermodynamic  equilibration  between 
species,  and  species  production  through  atomic  reactions,  e.g.  ionization,  recombina¬ 
tion,  charge  exchange.  The  Boltzmann  equation  coupled  with  Maxwell’s  equations  for 
electromagnetic  fields  completely  describe  the  plasma  dynamics.  Plasmas  have  been 
simulated  using  this  model  with  specific  forms  of  the  collision  operator  (e.g.  Vlasov 
equation  and  Fokker-Planck  equation).  [7-16]  However,  the  Boltzmann  equation  spans 
six  dimensions  corresponding  to  spatial  position  and  velocity,  in  addtion  to  time.  As 
a  consequence  of  the  large  dimensionality  plasmas  are  simulated  using  the  Boltzmann 
equation  only  when  required  to  capture  the  essential  physics.  The  applications  are 
generally  limited  to  plasmas  with  narrow  distributions,  small  spatial  extent,  and  short 
time  durations. 
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Particle  in  cell  (PIC)  plasma  models  apply  the  Boltzmann  equation  to  represen¬ 
tative  superparticles,  which  are  far  fewer  than  the  number  of  particles  in  the  actual 
plasma.  [17]  In  this  manner,  PIC  methods  provide  a  statistical  sampling  of  phase  space. 
While  PIC  methods  have  been  successful  in  modeling  many  physical  effects  [18-20], 
they  are  not  universally  applicable  due  to  grid  effects  and  statistical  errors  or  particle 
noise,  which  scales  with  the  number  of  particles  as  IW1/2  [21].  PIC  simulations  have 
similar  limitations  as  simulations  using  kinetic  theory. 

Another  approach  to  capture  more  complete  physics  is  to  generalize  the  single-fluid 
MHD  model  that  results  from  moments  of  the  Boltzmann  equation,  Eq.  (1).  The 
generalization  described  here  allows  for  multiple  species,  and  the  resulting  model  is 
the  multi-fluid  plasma  model.  Each  fluid  is  assumed  to  have  a  Maxwellian  velocity 
distribution.  The  generalization  also  allows  for  atomic  reactions  such  ionization,  re¬ 
combination,  and  charge  exchange.  Furthermore,  the  moment  equations  can  include 
higher  moments  to  more  accurately  model  the  evolution  of  plasmas  that  deviate  from 
thermodynamic  equilibrium. 

1.2  Introduction  to  Fluid  Plasma  Models 

Taking  moments  of  the  Boltzmann  equation,  Eq.  (1),  provides  equations  that  govern  the 
evolution  of  the  moment  variables.  The  moment  variables  are  defined  from  moments 
of  the  distribution  function.  For  example, 


na=  j  fa{v)dv, 

(2) 

nauai  =  f  vifa(v)dv , 

(3) 

where  the  integrals  are  performed  over  all  velocity  space  and  i  represents  the  spatial 
coordinate  index.  The  resulting  fluid  variables  are  number  density  na,  velocity  uQ.  ■, 
etc.  Moments  of  the  Boltzmann  equation  provide  evolution  equations  for  these  moment 
variables.  The  governing  equations  for  the  limiting  case  of  a  collisionless  plasma  with 
only  two  species,  ions  and  electrons,  the  two- fluid  plasma  model  is  presented  in  Ref.  [22]. 

The  governing  equations  of  the  two-fluid  plasma  model  can  be  combined  to  form 
the  single-fluid  MHD  model.  [23]  In  the  derivation  of  the  MHD  model  several  ap¬ 
proximations  are  made,  which  limit  its  applicability  to  low  frequency  phenomena  and 
ignores  potentially  significant  finite  electron  mass  and  charge  separation  effects.  These 
limitations  are  not  present  in  the  multi-fluid  plasma  model. 

Generalizing  the  moment  approach  to  include  an  arbitrary  number  of  species  and 
to  include  atomic  reactions  yields  the  multi-fluid  plasma  model.  The  derivation  follows 
that  presented,  for  example,  by  Braginskii  in  Ref.  [24],  However,  the  form  of  the 
equations  are  derived  here  for  the  conservation  variables  in  flux/source  form  where 
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hyperbolic  and  parabolic  fluxes  are  in  balance  with  source  terms.  The  equation  systems 
can  be  expressed  as 


d  ,JLF  -  o 

dtQa  +  dxk  ak  ~  a ’ 


(4) 


where  qa  is  the  vector  of  conservation  variables  of  species  a ,  FQ  is  the  tensor  of  hyper¬ 
bolic  fluxes  for  a,  and  Sa  is  the  source  vector  for  a.  Note  that  throughout  this  proposal 
vectors  and  tensors  are  represented  in  bold,  and  their  components  are  represented  in 
italics  with  subscript  indices  ( i,j,k,l ).  Repeated  indices  of  the  spatial  coordinate  are 
summed  in  the  usual  convention  of  Einstein  notation.  The  source  vector  includes  the 
coupling  to  the  other  fluid  species  and  to  the  electromagnetic  fields.  For  example,  the 
electric  and  magnetic  fields  appear  in  the  source  terms  of  any  charged-fluid  equations. 
The  field  dynamics  are  governed  by  Maxwell’s  equations,  which  have  source  terms  that 
contain  the  charged-fluid  variables.  See  Sec.  1.4.  Maxwell’s  equations  can  be  expressed 
in  the  form  given  by  Eq.  (4)  where  a  =  EM . 

Since  each  evolution  equation  derived  from  the  moment  approach  introduces  the 
next  higher  moment,  the  series  continues  indefinitely.  The  equation  system  must  be 
terminated  and  closure  relations  must  be  specified  that  relate  the  higher  moment  vari¬ 
ables  to  the  lower  moment  variables  in  the  system.  The  complete  multi-fluid  model 
and  its  extensions  are  presented  in  Sec.  1.3. 


1.3  The  Multi-Fluid  Plasma  Model 

The  governing  equations  for  multi-fluid  plasma  models  are  derived  by  taking  moments 
of  the  Boltzmann  equation,  Eq.  (1),  for  each  species,  as  briefly  introduced  in  Sec.  1.2. 
The  multi-fluid  plasma  model  (including  the  electromagnetic  equations)  are  expressed 
in  divergence  form  as  in  Eq.  (4).  Each  fluid  is  assumed  to  be  sufficiently  close  to 
thermodynamic  equilibrium  that  its  velocity  distribution  function  is  well  approximated 
by  a  limited  expansion  about  a  Maxwellian  distribution.  The  fluid  variables  are  derived 
by  taking  moments  of  the  distribution  function. 


Pa  =  TTla  j 

/a(v)dv, 

(5) 

Pa^ai  =  ?Tla  j 

[  Vifa(v)dv, 

(6) 

=  paTa  =  rna 

f  ^w2/«(w)dw, 

(7) 

where  ma  is  the  mass  of  species  a,  and  the  distribution  function  is  expressed  equiv¬ 
alently  as  a  function  of  either  the  velocity  v  or  the  random  velocity  about  the  mean 
fluid  velocity  w  =  v  —  uQ.  The  fluid  variables  for  each  species  are  mass  density  pa  (1 
component),  velocity  uQ  (3  components),  and  pressure  pa  (1  component).  Ta  is  the 
temperature  associated  with  species  a.  The  model  has  a  total  of  five  fluid  variables  or 
components  for  each  species  and  is  called  the  5 TV- moment  fluid  model. 
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An  evolution  equation  for  the  mass  density  of  each  species  is  given  by  the  zeroth 
moment  of  the  Boltzmann  equation. 


d_ 

at 


Pa  + 


d 

dxk 


(j PaUQk ) 


d 

mpa 


(8) 


The  net  mass  production  rate  of  species  a  due  to  atomic  reactions  is  denoted  on 
the  right-hand  side  of  the  equation  with  a  subscript  T.  Contributions  due  to  atomic 
reactions  are  described  later  in  Sec.  1.4.2. 

The  first  moment  of  the  Boltzmann  equation  yields  momentum  equations  and  de¬ 
scribes  the  evolution  of  the  momentum  density  for  each  species. 


iPa^-ai)  T  (Pa^ai  33ja^  T  Pa^ik)  QaTlct  (-Ej  T  ^ijk^ajRk)  ^  '  Rgfii  T  [Pa'U'a.i) 

(9) 

where  E  and  B  are  the  electric  and  magnetic  fields  and  Ra/j  is  the  momentum  transfer 
vector  from  species  a  to  species  (3  due  to  collisions.  I  is  the  identity  tensor.  The  number 
density  has  been  introduced  and  is  defined  by  na  =  pa/ma.  The  net  momentum 
production  rate  of  species  a  due  to  atomic  reactions  is  denoted  on  the  right-hand  side 
of  the  equation  with  a  subscript  T. 

The  second  moment  of  the  Boltzmann  equation  yields  an  energy  equation  for  each 
species,  which  is  expressed  in  divergence  form  for  the  total  energy. 


d  d 

“1“  qx  “1“  Pa)  ^ a^,  ^akl  =  Qa^a^aj Ej  + 


2_j  U°3 


Rafij  +  Qa/3  +  +  ~g^£o 


d 


(10) 

where  ha  is  the  heat  flux  vector,  Qap  is  the  heat  generated  in  species  a  due  to  collisions 
with  species  /3,  and  the  total  energy  is  defined  by 


£a  — 


7-1 


Pa  +  ^PaUQ. 


(11) 


and  7  is  the  ratio  of  specific  heats.  The  evolution  equation  for  the  total  energy  can 
be  combined  with  the  previous  two  moment  equations  to  provide  an  expression  for  the 
evolution  of  pressure.  The  energy  addition  rate  of  species  a  due  to  atomic  reactions  is 
denoted  on  the  right-hand  side  of  the  equation  with  a  subscript  T. 

Since  the  heat  flux  represents  a  higher  moment  of  the  distribution  function,  a  closure 
relation  must  be  specified  that  relates  it  to  the  lower  moment  variables.  Fourier’s  law 
is  a  commonly  used  relation,  which  gives  the  ith  component  of  the  heat  flux  as 


h  — 

,bai  —  o 

3  dxi 


(12) 


where  naij  is  the  ij  component  of  the  thermal  conductivity  tensor  which,  in  general, 
depends  on  the  strength  and  relative  orientation  of  the  magnetic  field.  Additional 
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closure  relations  are  needed  if  a  pressure  tensor  is  used  instead  of  the  scalar  pressure 
used  in  the  momentum  and  energy  evolution  equations. 

The  5  A-rnoment  model  and  appropriate  closures  are  derived  using  a  Chapman- 
Enskog  expansion  of  the  distribution  function.  The  distribution  function  is  assumed  to 
be  a  Maxwellian  (thermodynamic  equilibrium)  distribution  with  an  expansion  in  pow¬ 
ers  of  a  small  parameter  given  by  the  Knudsen  number,  the  ratio  of  the  mean  free  path 
to  the  characteristic  plasma  size.  Closure  relations  for  the  transport  coefficients  are 
found  by  retaining  only  the  linear  terms  of  the  expansion.  For  examples  of  calculations 
of  transport  coefficients,  see  Ref.  [24,25]. 


1.4  Maxwell’s  Equations 

The  electromagnetic  fields  influence  the  motion  of  the  plasma  fluid  through  the  Lorentz 
force,  which  is  contained  in  Eq.  (1)  for  the  kinetic  model,  Eq.  (9)  for  the  5  A-moment 
model,  and  Eq.  (43)  for  the  13 A- moment  model.  The  motion  of  the  plasma  influences 
the  evolution  of  the  electromagnetic  fields  through  the  redistribution  of  charge  density 
and  current  density.  Maxwell’s  equations  govern  the  evolution  of  the  electromagnetic 
fields.  The  net  charge  density  and  total  current  density  are  calculated  directly  from 
the  plasma  state,  i.e.  the  distribution  functions  or  the  multi-fluid  plasma  variables,  as 


Pc  =  ^2  qana  =  ^2  qa  /  fa  (v)g?v 

a  a 

(13) 

*  =  £&««««,  =£fe  [vifaMdv. 

a  a  ^ 

(14) 

These  terms  appear  as  source  terms  in  Maxwell’s  equations,  which  can  be  expressed  as 


d  d 

_  -eijkfaTEk 


dt 


d  d 

k  Qx  B k  MO  Ji 

0 

—  Pc 

OXk 

J-Bk  =  0 

axk 


(15) 

(16) 

(17) 

(18) 


The  divergence  constraint  relations  on  E  and  B  should  not  enter  the  calculation 
of  the  field  dynamics.  Mathematically,  if  the  initial  fields  satisfy  the  divergence  con¬ 
straints,  then  the  field  evolution  maintains  the  constraints.  Numerically,  the  divergence 
constraints  must  be  explicitly  enforced  by  “cleaning”  the  fields  with  either  an  elliptic 
method  [26],  a  hyperbolic  method  [27],  or  a  parabolic  method  [28]. 

The  Jacobians  of  the  hyperbolic  fluxes  dFa/dqa  of  the  governing  equations  are 
constructed  in  the  usual  way  from  Eq.  (4).  The  eigenvalues  of  the  flux  Jacobians  give 
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the  characteristic  velocities.  In  one  dimension,  the  eigenvalues  of  the  fluid  equations 
are 

A fluid  {vax !  vax  i  csa}  (19) 

where  the  acoustic  speed  for  species  a  is  defined  as 


The  electron  acoustic  speed  is  larger  than  the  ion  acoustic  speed  for  the  same  fluid 
temperatures  due  to  the  large  ion  to  electron  mass  ratio.  The  electron  acoustic  speed 
can  be  larger  than  the  Alfven  speed,  which  is  a  component  of  the  eigenvalues  of  MHD. 
The  Alfven  speed  for  an  ion/electron  plasma  is  defined  as 


where  B  is  the  magnitude  of  the  magnetic  field  and  is  the  permeability  of  free  space 
(47T  x  10“ 7  H/m).  The  eigenvalues  of  the  field  equations  are 

A  field  =  {±c}  ,  (22) 

where  c  is  the  speed  of  light.  Therefore,  the  eigenvalues  of  the  multi-fluid  plasma 
model  are  generally  not  bounded  by  the  eigenvalues  of  the  MHD  model.  In  general,  the 
fastest  times  that  must  be  resolved  in  the  multi-fluid  plasma  model  are  the  timescales 
associated  with  the  electromagnetic  fields  and  the  electron  fluid,  namely,  the  light 
transit  time  and  the  electron  plasma  oscillation  period. 

1.4.1  Collisional  Effects 

The  fluids  of  the  multi-fluid  plasma  model  interact  primarily  through  the  electromag¬ 
netic  fields,  which  produce  long-range  forces.  However,  short-range  collisional  effects 
can  have  a  significant  impact  on  the  overall  plasma  behavior  and  evolution.  Specifi¬ 
cally,  collisional  effects  can  thermalize  the  fluids  such  that  the  entire  plasma  approaches 
a  thermodynamic  equilibrium.  The  time  for  thermalization  within  a  species  is  the 
self-relaxation  time.  Collisions  between  species  bring  the  fluids  into  thermodynamic 
equilibration  with  a  characteristic  relaxation  time.  For  example,  an  electron-ion  fluid 
plasma  typically  has  relaxation  times  such  that  ree  <C  Tu  <C  Tje. 

Collisional  effects  are  formally  treated  in  the  fluid  model  by  evaluating  the  collision 
integral,  the  right-hand  side  of  Boltzmann  equation,  Eq.  (1).  The  Landau  form  of  the 
collision  integral  for  Coulomb  collisions  [29]  is  expressed  as 
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where 


Uki  = 


Vi  —  v' 


(Vj  -  v'j)2  6ki  -  (vk  -  v'k )  (vi  -  v[) 


(24) 


The  collisional  terms  accounts  for  transfer  of  momentum  and  energy  between  the 
different  fluids.  [24]  The  terms  appear  as  frictional  effects  in  Eqs.  (9)  and  (10),  Yhp  Rafii- 
In  general,  the  momentum  transfer  vector  is  defined  by  the  first  moment  of  the  collision 
integral  with  the  random  velocity,  as  used  in  Eq.(7), 


Rapi  —  vna 


UJiCapdw. 


(25) 


If  only  binary  collisions  that  result  in  small  angle  deflections  are  included,  the  rate  of 
momentum  transfer  from  species  a  to  species  (5  is  given  by 

m0 

Rah  =  - - - PaVap  ( Vai  -  Vfo)  ,  (26) 

ma  +  mp 

where  vap  is  the  collision  frequency  between  species  a  and  f3.  Momentum  conservation 
requires  Rapi  =  —RpQi  and  Raa,  =  0.  Heat  generation  due  to  random  collisions  is 
defined  by  the  second  moment  of  the  collision  integral  with  the  random  velocity, 

Qap  =  \m<*  J  w2Capdw 

These  relations  hold  for  nonreacting  collisions,  such  as  elastic  collisions.  Including 
atomic  reactions  requires  a  generalization  or  additional  effects  to  the  collisional  effects. 
Including  the  effect  of  atomic  reactions  on  the  plasma  dynamics  requires  additional 
terms  in  the  collision  operator,  and  is  described  in  the  next  section. 


(27) 


1.4.2  Atomic  Reactions 

The  purpose  of  including  the  effect  of  atomic  reactions  into  the  multi-fluid  plasma 
model  is  to  capture  the  time-dependent  ionization,  recombination,  and  charge  exchange 
reactions  that  are  important  in  laboratory  and  transient  plasmas.  Atomic  reactions 
lead  to  the  transfer  of  density,  momentum,  and  energy  between  the  fluids  in  the  multi¬ 
fluid  plasma  model.  For  example,  ionization  depletes  the  neutral  fluid  and  increases 
the  ion  and  electron  fluids. 

Contributions  from  atomic  reactions  are  identified  by  terms  on  the  right-hand  side 
of  Eqs.  (8),  (9),  and  (10).  The  terms  are  expressed  below  for  the  ionization  (ion), 
recombination  (rec),  and  charge  exchange  (cx)  reactions  for  a  hydrogen  plasma  com¬ 
posed  of  neutral  hydrogen  (a  =  n),  ionized  hydrogen  (a  =  i),  and  electrons  (a  =  e). 
Additional  reactions  and  other  plasma  constituents  are  also  possible. 

The  formal  procedure  to  account  for  atomic  reactions  in  the  fluid  model  requires 
evaluating  the  collision  integral,  the  right-hand  side  of  Boltzmann  equation.  The  eval¬ 
uation  involves  convolving  the  distribution  functions  of  the  reacting  species  with  the 
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reaction  cross  section  that  depends  on  their  relative  velocities.  For  example,  the  result¬ 
ing  effect  on  the  ion  species  from  charge  exchange  reactions  with  the  neutral  species 
is 

dfi  I 


dt 


=  /  °cx  |v  -  v|  [/j(v')/n(v)  -  /i(v)/n(v')]  dV, 


(28) 


where  acx  =  acx  (|v7  —  v|)  is  the  cross  section  for  charge  exchange.  [30] 

Evaluating  Eq.  (28)  is  accomplished  by  assuming  distribution  functions  for  each 
species,  for  example  a  Chapman-Enskog  expansion.  Zeroth  order  contributions  are 
found  by  assuming  a  Maxwellian  distribution.  These  are  presented  below.  More  de¬ 
tailed  calculations  will  be  performed  for  the  proposed  work  to  include  higher  order 
corrections  that  account  for  small  deviations  away  from  Maxwellian. 

The  contributions  to  the  species  densities  are 

d 


Pn 


dt 

d_ 

Ft 

d_ 

Ft 


—  ?nn  (  Tjon  T  Tr 


—  m i  (T ion  rrec)  , 


—  ITT-e  (P ion  rrec)  . 


The  contributions  to  the  species  momenta  are 

d 


dt 


( Pnuni , 


|  (««.) 


d_ 

Ft 


( PeUei ] 


—  (F ion  T  Pcx)  TTlnUni  -j-  (r rec  T  Fez)  Vlj lli^ , 


—  (r ion  T  r cx  )  n7  /  V./r//  (Free  T  Pea: )  , 


—  (Tjon  T  rca;)  meUni  (rrec  T  rca;)  1TleUe 


The  contributions  to  the  species  energies  are 

d 

Fte 

d_ 

Ft 


—  (r*on  t  rca;)  £n  +  (rrec  +  rca:)  £j, 


—  (Fjon  T  Fez)  (rrec  T  r^)  £j, 


d 

Ft£t 


—  (F*on  T  Fez)  £n  (Free  T  FCa;)  £e. 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


The  reaction  rates  are  given  by 

r ion  —  (vv)ion  nnrii 
Trec  =  ( av)recneni 
Fez  =  ( vv)cx  nnrii 

As  indicated  in  Eq.  (28),  the  cross  sections  a  and  reaction  rate  parameters  (c tv)  can 
depend  on  temperature. 


(38) 
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1.5  Higher  Moment  Models 

The  5  A-mornent  model  presented  above  provides  an  adequate  description  for  many 
plasmas  that  have  relaxed  sufficiently  close  to  thermodynamic  equilibrium.  For  truly 
non-equilibrium  plasmas,  a  kinetic  model  is  needed.  Generalizations  of  the  moment 
method  (higher  moment  models)  extend  the  applicability  of  fluid  models  to  non¬ 
equilibrium  plasmas.  The  higher  moments  extend  the  ability  to  characterize  the  dis¬ 
tribution  function  with  averaged  or  moment  variables. 

The  13 A- moment  model  provides  the  next  logical  extension  of  the  5 A- moment 
plasma  model.  The  model  as  first  introduced  by  Grad  [31]  and  has  since  been  refined  by 
many  others,  for  example  Refs.  [32-36].  Specifically,  Ref.  [32]  provides  a  derivation  for 
gas  dynamics  using  a  Pearson-type-IV  distribution  for  the  velocity  distribution  function 
that  better  matches  experimental  measurements  and  improves  the  hyperbolicity  of 
the  resulting  evolution  equations.  These  are  important  advantages  over  Grad’s  initial 
method.  The  results  are  only  briefly  presented  here  for  a  neutral  gas.  We  have  extend 
the  model  to  derive  a  13A-moment  fluid  model  for  a  multicomponent  plasma. 

As  before,  moments  of  the  distribution  function  provide  the  fundamental  fluid  vari¬ 
ables.  Consistent  with  the  original  derivation  by  Grad  [31],  only  the  variables  with  an 
explicit  physical  meaning  are  retained.  For  each  species  fluid  the  density  and  velocity 
is  defined  as  in  Eqs.  (5)  and  (6).  The  scalar  pressure  is  replaced  with  the  pressure 
tensor  and  the  heat  flux  vector  is  also  defined  as  a  solution  variable.  These  variables 
are  defined  as 


where  the  species  subscript  a  has  been  dropped  for  clarity.  In  addition  to  the  4  com¬ 
ponents  of  density  and  velocity,  the  model  has  a  total  of  thirteen  fluid  variables  or 
components  for  each  fluid.  The  pressure  tensor  p  has  6  components,  since  pt]  =  pp 
and  heat  flux  h  has  3  components.  Note  the  scalar  pressure  is  the  trace  of  the  pressure 
tensor,  and  is  given  by 


p  =  pT  =  m 


/(w)rfw, 


(41) 


which  is  consistent  with  the  definition  from  Eq.  (7). 

Taking  moments  of  the  Boltzmann  equation  provides  the  evolution  equations  for 
the  moment  variables. 


d  d  .  . 


d_ 

Ft 


( pui )  + 


d 

dxk 


{puiUk  +pik)  =  0, 


d_ 

Ft 


( pUiUj  +  Pi  j )  + 


d 

dxk 


( puiUjUk  +  3  pijUk  +  rriijk)  =  0, 


(42) 

(43) 

(44) 
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d_ 

dt 


Puj+Pjj 
2 

d 

dxk 


ii;  +  Pi ju j  +  hi 


+ 


pUj+Pjj 


UiUk  + 


PikV-i 


-(-  2 UiPjkUj  -|-  2h i>lk  TTl’ijkUj  -|- 


h* //.- 
~Y 


0.  (45) 


Equations  (44)  and  (45)  provide  evolution  equations  for  the  pressure  tensor  and  heat 
flux  vector.  The  equation  system  must  be  closed  for  the  variables  corresponding  to  the 
higher  moments,  namely, 


TTlijk 

J 

1  WiWjWkf(w)dw, 

(46) 

ii 

r  2 

\v~WiWjf(w)dw. 

(47) 

Closure  relations  are  derived  by  assuming  a  particular  form  of  the  distribution  function. 
For  the  proposed  work,  a  Pearson-type- IV  distribution  will  be  used,  as  in  Ref.  [32], 


1.6  Computational  Solution  Methods 

The  computational  methods  used  to  solve  the  multi-fluid  plasma  model  and  Maxwell’s 
equations  described  above  are  developed  to  be  optimally  matched  for  the  expected 
physics  for  each  set  of  governing  equations.  The  spatial  representation  for  all  of  the 
system  variables  is  based  on  finite  element  methods  where  the  simulation  domain  is 
divided  into  discrete  elements,  which  are  either  quadrilaterals  (2D)  or  hexahedrals  (3D). 
The  variation  of  the  solution  variables  within  each  element  are  modeled  by  projecting 
the  variables  onto  a  set  of  spatially  dependent  basis  functions  Vh  of  order  h,  such  that 
within  each  element  Q  variable  q  is  represented  as 

9n(x)  =  ^2qnhvh(x)-  (48) 

h 

The  basis  functions  are  Legendre  polynomials  or  Jacobi  polynomials  for  the  proposed 
implementation.  The  finite  element  representation  captures  high-order  spatial  varia¬ 
tions,  which  is  critical  when  anisotropic  properties  exist  [37]  or  when  hyperbolic  fluxes 
are  balanced  by  the  fluxes  from  other  species,  which  occurs  at  equilibrium. 

The  governing  equations  for  the  multi-fluid  plasma  model  can  be  expressed  as  given 
by  Eq.  (4).  The  equations  can  be  grouped  according  to  their  physically  characteris¬ 
tic  spatial  and  temporal  scales,  which  forms  the  basis  for  selecting  the  appropriate 
computational  method.  The  fast  temporal  and  short  spatial  dynamics  of  the  electron 
fluid  and  electromagnetic  fields  suggests  solving  these  governing  equations  in  a  coupled 
manner  as 


d 

qem 

9 

F 'EMk 

Sem 

dt 

.  . 

+  dxk 

.  Fek  . 

.  Se  . 
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where  the  source  terms,  in  general,  depend  on  the  solution  variables  of  the  fields  and 
other  species,  Sa  =  Sa  ( qEM,Qe,Qi,Qn )•  The  slow  temporal  and  long  spatial  dynamics 
of  the  ion  and  neutral  fluids  suggests  solving  these  governing  equations  in  a  separately 
coupled  manner  as 


di 

d 

-\ - 

'si 

Qn_ 

dxk 

J^nk. 

_Sn_ 

(50) 


where  the  source  terms  again,  in  general,  depend  on  the  solution  variables  of  the  fields 
and  other  species. 

A  Galerkin  method  is  used  to  obtain  spatially  discretized  equations.  The  governing 
equations  are  multiplied  by  each  basis  function  and  integrated  over  the  element  volume 
to  give  the  integral  equation 


Lv"BidV  + 


an 


vhFkdAk  -  [  Fk-^—vhdV  =  [  vhSdV. 
Jn  oxk  Jn 


(51) 


The  integrals  are  evaluated  by  Gaussian  quadrature.  The  source  terms  are  also  pro¬ 
jected  onto  the  basis  functions,  similar  to  Eq.  (48).  The  solution  of  Eq.  (51)  throughout 
the  spatial  domain  and  its  time  advance  provides  the  complete  time-dependent  solu¬ 
tion.  Both  the  spatial  representation  and  time  advance  are  optimized  based  on  the 
expected  physics. 


1.7  Implementation  and  Computational  Demonstrations 

This  section  provides  a  brief  description  of  the  implementation  details  and  presents 
demonstrations  of  the  computational  algorithm  generated  by  our  Computational  Plasma 
Dynamics  Group  at  the  University  of  Washington.  We  have  developed  accurate  numer¬ 
ical  algorithms  for  the  two-fluid  plasma  model  that  are  based  on  finite  volume  [22,38] 
and  DG  methods  [39-42].  Higher  moment  plasma  models  (10  A  and  13 N)  have  been 
derived  and  investigated.  [34-36,  43]  Numerical  formulations  have  been  investigated 
that  allow  for  the  implicit  solutions  of  plasma  models.  [2]  These  computational  ad¬ 
vances  have  been  applied  to  study  plasma  phenomena  with  high  fidelity  models,  e.g. 
Refs.  [44-46].  In  addition,  we  have  developed  non-reflecting  open  boundary  conditions 
that  allow  simulating  infinite  space  on  finite  domains.  [47, 48]  The  results  described 
here  focus  on  our  developments  of  the  FE  method  that  we  have  applied  to  contin¬ 
uum  plasma  models  and  on  our  implementation  of  appropriate  boundary  conditions. 
The  algorithm  research  and  development  have  been  implemented  in  a  flexible  soft¬ 
ware  framework  called  WARPX  (Washington  Approximate  Riemann  Plasma),  which 
uses  C++  object  oriented  programming  and  other  modern  software  techniques  to  sim¬ 
plify  the  maintainability  and  extensibility  of  the  code  and  HDF5  for  parallel  output. 
WARPX  uses  MPI  message  passing  for  parallel  computer  architectures  and  OpenCL 
for  GPU/many-core  computer  architectures. 


13 


Physics-Based  Multi-Fluid  Plasma  Algorithm 


Shumlak 


1.7.1  Multi-Fluid  Plasma  Model 

We  have  developed  a  DG  method  [49-51]  to  solve  the  governing  equations  of  the  multi¬ 
fluid  plasma  model  on  a  computational  grid.  The  DG  method  is  a  finite  element 
approach  that  allows  for  arbitrarily  high-order  basis  functions  to  model  the  variation 
of  the  system  variables,  as  described  in  Sec.  1.6  and  in  more  detail  in  the  references 
provided. 

The  surface  integral  in  Eq.  (51)  uses  numerical  hyperbolic  fluxes  computed  with 
a  Roe-type  approximate  Riemann  solver  [22,39,52]  or  with  a  Lax-Friedrich  flux  [40]. 
The  overall  solution  is  built  upon  the  solutions  to  the  Riemann  problem  defined  by 
the  discontinuous  jumps  in  the  solution  at  each  element  interface.  Continuity  of  the 
fluxes  between  the  elements  ensures  conservation.  The  numerical  flux  for  a  first-order 
accurate  (in  space)  Roe-type  solver  is  written  in  symmetric  form  as 

P).+ 1/2  =  2  (-^>+1  A  -^i)  —  7 j  ^  Ik  (tfi+1  —  Qi )  |^fc|  (52) 

k 

where  r*,  is  the  kth  right  eigenvector,  is  the  kth  eigenvalue,  and  A  is  the  kth  left 
eigenvector,  evaluated  at  the  element  interface  (i  +  1/2).  The  values  at  the  element 
interface  are  obtained  by  a  Roe  average  of  the  neighboring  elements.  The  flux  calcu¬ 
lated  as  above  is  normal  to  the  element  interface  which  is  the  desired  orientation  for 
calculating  the  surface  integral. 

As  a  simple  example,  a  the  two-dimensional,  second-order  accurate  algorithm,  uses 
a  set  of  linear  basis  functions. 


{vh}  =  {vo,vx,vy} 


(  x-Xij  y  -  yij  \ 
I  ’  Ax/2  ’  Ay/2  j 


(53) 


where  the  center  of  the  mesh  element  is  located  at  ( Xij,yij )  and  extends  Ax  by  Ay. 
The  conserved  variables  q  are  defined  as 


q  =  qo  +  qxvx  +  qyvy  (54) 

within  each  mesh  element.  The  DG  method  has  been  implemented  using  a  modal 
approach,  as  opposed  to  the  nodal  approach  that  is  common  in  typical  FE  methods. 
Update  equations  for  the  coefficients  for  each  conserved  variable  are  found  directly 
from  Eq.  (51)  applied  to  each  mesh  element.  The  temporal  evolution  is  computed  with 
a  Runge-Kutta  method.  A  third  order  TYD  method  has  been  used  successfully.  [53] 
High-order  FE  methods  more  accurately  capture  fine  scale  structures,  which  is 
demonstrated  by  studying  the  propagation  of  dispersive  waves.  Electrostatic  ion  cy¬ 
clotron  waves  produce  dispersive  waves  which  have  an  analytical  description.  A  wave 
with  frequency  uj  and  wave  number  k  is  described  by  the  dispersion  relation,  u2  = 
k2c2  +  lo2,  where  cs  is  the  sound  speed  and  luc  is  the  cyclotron  frequency.  Nine  modes 
are  excited,  such  that  the  analytic  solution  evolves  to  contain  many  frequencies  and 
wave  numbers,  as  seen  in  Fig.  1.  We  simulated  the  problem  using  the  DG  method 
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Figure  1:  Dispersive  electrostatic  ion  cyclotron  waves  generated  from  a  nine  mode  excitation. 
Plotted  is  the  analytical  solution  and  numerical  solutions  with  a  fixed  number  of  unknowns: 
3rd  order  DG  with  33  elements  (dg-3),  8th  order  DG  with  12  elements  (dg-8),  and  16th  order 
DG  with  6  elements  (dg-16).  The  benefit  of  high-order  is  evident. 


with  different  orders  of  accuracy  but  with  the  same  effective  resolution  (approximately 
constant  number  of  unknowns).  [42]  The  results  in  Fig.  1  show  that  a  high-order  DG 
method  better  matches  the  analytical  solution. 

High-order  FE  methods  also  better  resolve  anisotropic  behavior  that  often  results 
in  magnetized  plasmas.  The  ratio  of  parallel  to  perpendicular  thermal  conductivities 
can  be  as  large  as  106.  Low-order  methods  introduce  significant  numerical  diffusion 
that  obscures  the  anisotropic  transport  properties.  We  have  investigated  the  ability  of 
high-order  FE  to  resolve  perfect  heat  conduction  along  a  toroidal  magnetic  field  in  a 
Cartesian  grid,  where  the  ratio  of  the  conductivities  is  D^_/D^  =  0.  [37]  The  results  are 
presented  in  Fig.  2,  which  show  that  higher-order  FE  maintains  a  lower  numerically 
measured  D±  even  for  the  same  effective  resolution. 

The  discontinuous  Galerkin  algorithm  has  been  applied  to  the  electromagnetic 
plasma  shock  demonstrating  the  transition  from  gas  dynamic  shocks  to  the  MHD 
shock  [54, 55]  as  the  Larrnor  radius  is  reduced.  Analysis  of  the  data  shows  the  differ¬ 
ences  caused  by  the  additional  plasma  waves  that  are  captured  in  the  two-fluid  model 
and,  consequently.  [22]  It  also  illustrates  the  dispersive  nature  of  the  waves,  which 
makes  capturing  the  effect  difficult  in  MHD  algorithms.  The  electromagnetic  plasma 
shock  serves  to  validate  the  algorithm  to  published  data  (MHD  limit)  and  analytical 
results  (gas  dynamic  limit).  The  algorithm  has  also  been  applied  to  study  collisionless 
reconnection  and  the  results  are  compared  to  published  results  of  the  GEM  challenge 
problem.  [56]  The  problem  is  difficult  to  model  and  provides  a  rigorous  test  for  the  algo¬ 
rithm  and  benchmarks  to  other  algorithms.  The  evolution  of  the  reconnected  magnetic 
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Figure  2:  Numerically  measured  thermal  conductivity  perpendicular  to  a  toroidal  magnetic 
field.  Analytically,  D±  should  be  zero.  Plotted  are  the  values  of  D±  for  different  element  size 
(h)  and  FE  order  (np).  Higher-order  FE  maintains  a  lower  D  j_  even  for  the  same  effective 
resolution  ( dof ).  The  benefit  of  high-order  is  evident. 


flux  compares  remarkably  well  with  the  published  data.  [39] 

Additional  applications  of  the  5  A-moment  plasma  fluid  model  have  investigated 
lower  hybrid  drift  instabilities  to  explain  anomalous  resistivity  observed  in  experimental 
plasmas  [45]  and  the  evolution  of  FRC  plasmas  [44].  Applications  have  also  isolated 
and  contrasted  the  physics  captured  by  the  multi-fluid  plasma  model  and  Hall-MHD 
model.  [46] 

We  have  also  developed  advanced  physics  models  that  describe  plasma-neutral 
interactions  in  a  computationally  tractable  manner.  [57]  The  addition  of  a  neutral 
fluid  has  enabled  modeling  of  plasma  interactions  with  neutral  flows  and  sheath  for¬ 
mations.  [41,58]  These  simulations  have  demonstrated  the  generation  of  a  nonlinear 
Langmuir  wave  that  propagates  away  from  the  electrodes. 

The  13  A-moment  plasma  fluid  model  represents  a  unique  advance  because  it  extends 
the  region  of  validity  of  moment  models  towards  the  kinetic  regime.  In  this  manner, 
it  provides  a  meaningful  step  towards  bridging  the  gap.  We  have  developed  forms  on 
the  13 A- moment  plasma  fluid  model  based  on  the  Pearson- type- IV  distribution  for  the 
velocity  distribution  function.  Our  initial  results  in  Fig.  3  illustrate  the  relationship 
between  the  5A-moment  model,  the  13A-moment  model,  and  the  ooA-moment  model 
(continuum  kinetic  model) .  The  kinetic  solution  is  the  numerical  solution  of  the  Boltz¬ 
mann  equation  with  a  BGK  collision  operator.  A  similar  collision  operator  is  required 
for  the  13A-moment  model  to  better  approximate  the  physical  processes. 
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Figure  3:  Shock  tube  simulation  with  three  different  plasma  models:  5 A-morrient  model, 
13 A- moment  model,  and  ocA-moment  model  (continuum  kinetic  model).  The  comparison 
illustrates  the  convergence  and  relationships  between  the  models. 
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Figure  4:  Nonlinear  solution  of  the  Vlasov-Poisson  equation  system  for  the  strong  Landau 
damping  problem  using  seventh-order  polynomials  on  a  20  x  80  phase-space  grid  in  (x,vx) 
at  t  —  60 Up1.  Evolution  of  the  potential  energy  exhibits  the  expected  damping  and  growth 
behavior. 

1.7.2  Continuum  Kinetic  Plasma  Model 

As  indicated  by  the  oo-Moment  solution  in  Fig.  3,  we  have  made  substantial  progress 
in  solving  continuum  kinetic  plasma  models  -  directly  evolving  Boltzmann  equation, 

Eq.  (1).  Using  the  DG  method  to  solve  the  Vlasov-Poisson  system  has  recently  gained 
substantial  interest.  [59,  60]  We  have  also  implemented  the  DG  method  to  solve  the 
Vlasov-Poisson  equation  system  in  WARPX.  Figure  4  shows  the  nonlinear  solution  for 
the  strong  Landau  damping  problem  using  seventh-order  polynomials.  The  potential 
energy  exhibits  the  expected  damping  and  growth  behavior,  and  agrees  with  previously 
published  results.  [59,  61]  These  are  preliminary  calculations,  but  they  illustrate  the 
broader  applicability  of  the  DG  method.  In  particular,  high-order  representation  is  able 
to  accurately  capture  the  striations  that  occur  in  phase  space.  We  have  investigated  the 
conservation  properties  of  our  DG  method  implementation,  and  mass,  momentum,  and 
energy  are  conserved  significantly  better  than  in  other  continuum  implementations. 

In  addition,  we  have  also  developed  benchmark  problems  for  continuum  kinetic 
modeling.  While  several  benchmark  problems  exist,  they  are  limited  to  electrostatic 
and  two  dimensions.  [59,  61]  We  have  developed  a  three  dimensional  benchmark  for 
magnetized  plasmas.  [62]  It  has  been  initially  implemented  using  a  high-order,  finite- 
volume  method,  but  it  will  also  be  used  for  our  DG  implementation. 

1.7.3  Implicit  Time  Advance  and  the  Blended  Finite  Element  Method 

As  described  in  Secs.  1.6  and  1.7.1,  both  CG  and  DG  hnite  element  methods  have 
been  successfully  used  to  solve  continuum  plasma  models.  The  numerical  solution  of 
the  governing  equations  has  produced  some  unique  and  difficult  challenges.  Specifically, 
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the  advanced  plasma  models  encompass  a  wide  disparity  of  spatial  and  temporal  scales 
that  are  fundamentally  generated  by  the  massive  ions  and  neutrals,  lighter  electrons, 
and  massless  electromagnetic  fields.  Solving  the  multi-fluid  model  requires  a  numerical 
algorithm  that  can  capture  these  multiple  scales  efficiently  and  accurately. 

Building  on  the  success  of  the  explicit  time  advance  methods,  it  is  logical  to  inves¬ 
tigate  implicit  time  advance  methods  to  alleviate  the  time-step  limitations.  Implicit 
methods  allow  for  time  steps  that  are  not  limited  by  the  fast  waves.  The  entire  system 
of  governing  equations,  Eq.  (4),  is  rewritten  as 


for  the  electromagnetic  fields,  electron  fluid,  and  ion  fluid,  where  the  right  hand  side  of 
Eq.  (55)  includes  the  hyperbolic  fluxes  and  sources.  The  time  advance  can  be  expressed 
in  an  implicit  form  of  arbitrary  order  accuracy.  A  second-order  Crank-Nicolson  method 
has  been  implemented,  where  the  time  advance  equation  is  written  as 

qn+1=qn  +  ^[f(qn+1)  +  f(qn)],  (56) 

which  is  solved  iteratively  using  a  variety  of  implicit  numerical  methods. 

The  5 TV-moment  multi-fluid  plasma  model  contains  18  variables  for  the  ion  and 
electron  fluids  and  the  electromagnetic  fields.  Each  additional  fluid  requires  an  addi¬ 
tional  5  variables.  Formulating  an  implicit  time-advance  method  results  in  a  large,  stiff 
equation  system  with  a  matrix  that  is  difficult  to  invert.  We  tested  a  variety  of  im¬ 
plicit  numerical  methods,  e.g.  Conjugate  Residual,  Biconjugate  Gradient,  and  GMRES 
with  preconditioners  such  as  Incomplete  LU  and  Additive  Schwarz.  None  have  per¬ 
formed  adequately  without  artificially  altering  the  physics  of  the  problem  with  either 
unrealistically  massive  electrons,  e.g.  mi/me  =  25,  or  unphysically  slow  light  speed, 
e.g.  c/va  =  10.  The  difficulty  stems  from  the  widely  different  masses  of  the  ions  and 
neutrals  compared  to  the  electrons  and  fields.  For  example,  adequately  resolving  the 
electrons  in  space  and  time  results  in  greatly  “over-resolving”  the  ions  and  leads  to 
numerical  dissipation  of  the  ion  features. 

This  motivated  our  physics-based  computational  algorithm,  which  has  been  tremen¬ 
dously  successful.  The  algorithm  performs  a  physics-based  splitting  of  the  governing 
equations  and  applies  appropriate  spatial  and  temporal  representations  the  split  system 
while  still  providing  the  necessary  coupling. 

Because  of  the  light  mass  and  high  mobility  of  electrons,  electron  characteristic 
speeds  (eigenvalues)  are  fast  compared  to  ions  and  neutrals.  Electrons  only  form  true 
shocks  in  extremely  rare  situations,  for  example,  when  the  electron  or  ion  flow  veloc¬ 
ities  exceed  the  electron  thermal  speed,  which  does  not  usually  occur  since  the  ions 
generally  move  slower  than  electrons  and  the  ions  slow  the  electron  motion.  Ions  and 
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Figure  5:  Plasma  shock  simulation  using  the  two-fluid  plasma  model  with  physical  values 
for  me,  rrii,  and  c.  The  slow  compound  wave  structure  (in  blue  circle),  which  is  composed  of 
a  shock  and  rarefaction  wave,  produces  a  shock  in  the  ion  fluid  but  not  the  electron  fluid. 
Normalized  number  density  is  plotted  as  a  function  of  position. 

neutrals  have  much  slower  characteristic  speeds,  and  they  can  and  do  form  shocks.  The 
electromagnetic  coupling  between  the  ions  and  electrons  generates  forces  that  modify 
the  electron  spatial  structure  to  approximate  the  ion  spatial  structure.  If  the  ions  form 
a  shock,  the  electron  spatial  structure  has  a  variation,  but  typically  the  variation  is  not 
a  true  discontinuity.  An  example  of  this  behavior  is  shown  in  Fig.  5,  which  shows  a 
plasma  shock  simulation  using  the  two-fluid  plasma  model.  A  “slow  compound”  (SC) 
wave  structure  (x  =  0.23)  has  formed  at  the  base  of  the  rarefaction.  The  SC  wave 
structure  is  composed  of  a  shock  and  rarefaction  wave.  The  SC  structure  produces  a 
shock  in  the  ion  fluid,  but  the  electron  fluid  has  only  a  smooth  variation.  Since  the 
electromagnetic  fields  have  the  fastest  characteristic  speeds  of  the  system,  the  fields  in 
the  multi-fluid  plasma  model  should  never  form  shocks. 

We  have  developed  a  blended  finite  element  (BFE)  method  that  applies  a  high- 
order,  continuous  spatial  representation  (CG  method)  for  the  electron  fluid  and  elec¬ 
tromagnetic  fields  and  a  high-order,  discontinuous  spatial  representation  (DG  method) 
for  the  ion  and  neutral  fluids.  The  electron  fluid  and  electromagnetic  fields  do  not  re¬ 
quire  inter-element  flux  calculations  or  limiting  since  their  solutions  are  expected  to  be 
smooth  and  well-resolved.  The  ion  and  neutral  fluids  use  approximate  Riemann  fluxes 
and  limiting  to  capture  shocks.  Time  advance  of  the  electron  fluid  and  electromag¬ 
netic  fields  can  be  solved  implicitly  so  that  their  dynamics  do  not  limit  the  time  step. 

Since  the  dynamics  of  the  ion  and  neutral  fluids  typically  set  the  timescale  of  interest, 
their  evolution  is  solved  explicitly  using  the  TVD  Runge-Kutta  method  mentioned  in 
Sec.  1.6. 

The  performance  of  the  BFE  method  is  demonstrated  by  applying  it  to  study  laser- 
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Figure  6:  Blended  finite  element  method  applied  to  the  species  separation  problem  in  capsule 
implosions.  Number  densities  and  electric  field  are  shown  after  the  laser  drive  has  compressed 
the  multi-fluid  plasma  and  caused  a  species  separation.  Excellent  agreement  is  found  between 
the  explicit  and  implicit  solutions,  CFL=1  (solid  lines)  and  CFL=10  (dotted  lines  displaced 
vertically  for  clarity). 
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driven  implosions  to  compress  deuterium-tritium  fuel  to  fusion  conditions.  Separation 
of  the  deuterium  and  tritium  species  could  explain  the  experimentally  observed  low 
fusion  yield.  The  simulation  is  performed  using  the  BFE  method  with  the  multi-fluid 
5 TV- moment  plasma  model  for  three  fluids  (deuterium,  tritium,  and  electron).  The 
number  densities  are  shown  in  Fig.  6  after  the  laser  drive  has  started  the  compression. 
A  separation  clearly  develops.  The  solution  is  found  using  an  explicit  advance  (CFL=1) 
for  the  entire  solution,  qem ,  qe,  Qd,  Qt-  The  solution  is  also  found  using  an  implicit 
advance  (CFL=10)  for  the  electron  fluid  and  electromagnetic  fields,  qEMiQe,  and  an 
explicit  advance  for  the  ionic  species,  qD,QT-  The  agreement  can  be  seen  in  Fig.  6. 


1.7.4  Divergence  Cleaning  in  Maxwell’s  Equations 

Section  1.4  presented  Maxwell’s  equations.  The  time-dependent  equations,  Eqs.  (15) 
&  (16),  completely  describe  the  evolution  of  the  electromagnetic  fields.  Provided  the 
divergence  equations,  Eqs.  (17)  &  (18),  are  initially  satisfied,  the  evolved  fields  will 
always  satisfy  them  analytically.  However,  computational  solutions  can  lead  to  viola¬ 
tions  of  the  divergence  constraints  (involutions)  and  corresponding  nonphysical  effects 
such  as  charge  generation  and  parallel  magnetic  forces. 

We  have  previously  investigated  several  methods  to  “clean”  the  divergence  errors. 
These  methods  have  included  solving  the  mixed  potential  formulation  of  Maxwell’s 
equations  and  modifying  Maxwell’s  equations  to  be  purely  hyperbolic  [27].  Recently, 
we  have  investigated  a  parabolic  modification  [28]  that  produces  significantly  improved 
results  without  the  computational  expense  of  a  mixed  potential  formulation. 

In  this  approach,  the  divergence  constraint  relations  are  incorporated  with  the  time- 
dependent  field  equations  as  parabolic  terms.  Faraday’s  law  is  transformed  by  adding 
a  diffusive  term  proportional  to  the  divergence  error,  which  is  written  in  vector  form 
as 

HD 

_  =  _V  x  E  +  XBV  (V  •  B)  (57) 

where  xb  is  a  constant  magnetic  held  divergence  error  diffusivity.  The  evolution  of  the 
divergence  error  of  the  magnetic  held  is  found  by  taking  the  divergence  of  Eq.  (57)  to 
give 

^  (V  •  B)  =  XBV2  (V  •  B) ,  (58) 

which  locally  diffuses  the  divergence  error.  A  similar  modification  of  Ampere’s  law  is 
performed  to  give 


<9E 


eo/h)  -7^r  =  VxB-/ioXM“  +  e°V  ‘  E  -  XI  qaUa 


(59) 


where  xe  is  a  constant  electric  held  divergence  error  diffusivity. 

The  resulting  parabolic  equations  reduce  the  divergence  errors  throughout  the  do¬ 
main  at  any  given  moment.  Specifically,  the  error  is  not  moved  from  one  region  to 
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Figure  7:  Divergence  cleaning  comparison  between  the  parabolic  (blue)  and  hyperbolic  (red) 
methods.  Left  column  shows  the  solutions  at  t  =  0,  12,  24,  120.  Right  column  shows  the 
spatial  frequency  spectrum  of  the  divergence  error. 

another  as  it  is  with  the  hyperbolic  formulation.  The  diffusivities,  \B  and  xe,  can  be 
set  to  convenient  values  to  reduce  the  divergence  errors  to  acceptable  levels  without 
overly  restricting  the  time  advance.  Implicit  formulations  are  also  possible. 

The  parabolic  cleaning  method  has  been  implemented  in  ID  to  demonstrate  its 
improved  performance  compared  to  the  hyperbolic  formulation.  Figure  7  shows  a 
comparison  to  the  hyperbolic  cleaning  method.  The  spatial  frequency  content  of  the 
error  demonstrates  the  local  reduction  of  the  parabolic  cleaning  method.  The  method 
monotonically  removes  the  divergence  error  at  all  frequencies.  The  hyperbolic  method 
leads  to  a  broadband  spectrum  with  less  error  reduction. 

1.7.5  Open  Boundary  Conditions 

The  continuum  kinetic  and  multi-fluid  plasma  models  have  waves  with  disparate  speeds 
-  from  fast  light  waves  to  slow  ion  acoustic  waves.  Plasma  simulations  typically  have 
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timescales  of  interest  that  are  much  longer  than  the  transit  times  of  the  fastest  waves 
across  the  computational  domain.  If  the  simulation  is  modeling  an  unbounded  re¬ 
gion,  these  wave  should  leave  the  computational  domain  without  reflection.  However, 
many  numerical  treatments  for  open  boundary  conditions  generate  non-physical  re¬ 
flections  which  contaminate  the  solution.  We  have  implemented  non-local  boundary 
conditions  that  use  a  lacuna-based  method  [63]  where  an  overlapping  auxiliary  domain 
is  appended  to  the  interior  domain  and  the  boundary  condition  is  replaced  with  an 
interface  condition. 

The  governing  equation  in  the  interior  domain  is  expressed  as 

!9+a!Ft(<l)  =  s(‘l)'  (60) 

like  Eq.  (4),  where  a  has  been  dropped  for  clarity.  In  the  auxiliary  domain  the  governing 
equation  is  defined  as 

a  a 

+  g^Fk(w)  =  S(w) +  Q(q),  (61) 

where  £l(q)  is  the  “near-boundary  source”.  The  solutions  are  matched  in  a  transition 
region,  so  w  =  The  smooth  function  /i(x)  is  zero  at  the  beginning  of  the 

auxiliary  domain  and  one  at  the  interior /exterior  interface.  The  boundary  condition 
for  the  interior  solution  is  then  set  such  that  the  interface  condition  is 


<7 


=  W 

interface 


+ 

interface 


(62) 


The  auxiliary  solution  is  periodically  re-integrated  to  damp  the  leading  edge  of  the 
solution  before  it  reflects  and  contaminates  the  interior  solution. 

We  have  successfully  implemented  advanced  methods  for  numerically  treating  open 
boundaries  in  Refs.  [47,48].  The  lacuna-based  method  works  for  oblique  incidence 
waves  in  either  purely  hyperbolic  or  mixed  hyperbolic/parabolic  systems.  The  method 
even  works  in  two  dimensions  where  there  is  no  true  lacuna.  (Huygens’  principle  states 
true  lacunae  only  exist  in  odd  dimensional  space.)  Sample  results  are  shown  in  Fig.  8 
where  a  pressure  pulse  interacts  with  open  boundaries  on  all  sides  with  a  high  extinction 
rate  of  the  reflected  wave. 


2  Project  Personnel 

This  project  was  performed  by  Prof.  Uri  Shumlak  and  graduate  students  Eric  Meier, 
Andrew  Ho,  Robert  Lilly,  Sean  Miller,  Noah  Reddell,  Eder  Sousa,  and  Bhuvana  Srini- 
vasan.  Archival  journal  and  conference  papers  were  published  reporting  on  the  work 
from  this  project: 

•  G.V.  Vogman,  P.  Colella,  and  U.  Shumlak.  “Dory-Guest-Harris  instability  as 
a  benchmark  for  continuum  kinetic  Vlasov-Poisson  simulations  of  magnetized 
plasmas.”  Journal  of  Computational  Physics  277,  101-120  (2014). 
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Figure  8:  Lacuna-based  open  boundaries  interact  with  a  cylindrically  expanding  pressure 
pulse  and  exhibits  a  high  extinction  rate  of  the  reflected  wave. 

•  E.  Kansa,  U.  Shumlak,  and  S.  Tsynkov.  “Discrete  Calderon’s  projections  on 
parallelepipeds  and  their  application  to  computing  exterior  magnetic  fields  for 
FRC  plasmas.”  Journal  of  Computational  Physics  234,  172-198  (2013). 

•  E.T.  Meier  and  U.  Shumlak.  “A  general  nonlinear  fluid  model  for  reacting  plasma- 
neutral  mixtures.”  Physics  of  Plasmas  19(7),  072508  (2012). 

•  E.T.  Meier,  A.H.  Glasser,  V.S.  Lukin,  and  U.  Shumlak.  “Modeling  open  bound¬ 
aries  in  dissipative  MHD  simulations.”  Journal  of  Computational  Physics  231(7), 

2963-2976  (2012) 

•  U.  Shumlak,  R.  Lilly,  N.  Reddell,  E.  Sousa,  and  B.  Srinivasan,  “Advanced  physics 
calculations  using  a  multi-fluid  plasma  model.”  Computer  Physics  Communica¬ 
tions  182,  1767-1770  (2011). 

•  B.  Srinivasan  and  U.  Shumlak.  “Analytical  and  computational  study  of  the  ideal 

full  two-fluid  plasma  model  and  asymptotic  approximations  for  Hall-magnetohydrodynamics.” 
Physics  of  Plasmas  18(9),  092113  (2011). 

•  B.  Srinivasan,  A.  Hakim,  and  U.  Shumlak,  “Numerical  methods  for  two-fluid 
dispersive  fast  MHD  phenomena.”  Communications  in  Computational  Physics 
10(3),  183-215  (2011). 

•  J.  Loverich,  A.  Hakim,  and  U.  Shumlak,  “A  discontinuous  Galerkin  method  for 
ideal  two-fluid  plasma  equations.”  Communications  in  Computational  Physics  9, 

240-268  (2011). 

Dissertations  and  theses  can  be  obtained  from  the  University  of  Washington  library 
system  or  from  the  project  website,  http:/ /www. aa.washington.edu/research/cpdlab/. 
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Archival  publications  can  be  downloaded  from  ResearchGate. 

3  Conclusions 

Investigating  advanced  plasma  models  is  motivated  by  the  need  to  simulate  complicated 
plasma  physics  phenomena  that  is  not  captured  in  simpler  models.  The  multi-fluid 
plasma  model  is  proving  to  be  a  model  that  is  significantly  more  advanced  and  com¬ 
plete  than  the  usual  MHD  model.  Developing  numerical  algorithms  to  solve  these  more 
complete  plasma  models  that  exploit  the  different  physics  offers  computationally  effi¬ 
cient  methods  that  maintain  high  order  spatial  accuracy.  The  methods  extend  to  other 
continuum  desriptions,  such  as  continuum  kinetic  plasma  model  and  single-fluid  MHD 
models.  The  algorithm  developed  in  the  project,  its  implementation  into  WARPX,  and 
its  application  to  benchmark  and  real  experimental  problems  have  demonstrated  the 
capability  of  both  the  multi-fluid  plasma  model  and  the  numerical  techniques  in  the 
algorithm. 
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